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APPENDIX  D.  DERIVATION  OF  AND  02  INTEGRALS 


ABSTRACT 


Tracking  performance  depends  upon  the  quality  of  the 
measurement  data.  In  the  Kalman-Bucy  filter  and  other  trackers, 
this  dependence  is  well-understood  in  terms  of  the  measurement 
noise  covariance  matrix,  which  specifies  the  uncertainty  in  the 
values  of  the  measurement  inputs.  The  measurement  noise  and 
process  noise  covariances  determine  via  the  Riccati  equation,  the 
state  estimation  error  covariance.  When  the  origin  of  the 
measurements  is  also  uncertain,  one  has  the  widely-studied 
problem  of  data  association  (or  data  correlation) ,  and  tracking 
performance  depends  critically  on  additional  parameters, 
primarily  the  probabilities  of  detection  and  false  alarm.  In 
this  paper  we  derive  a  modified  Riccati  equation  that  quantifies 
(approximately)  the  dependence  of  the  state  error  covariance  on 
these  parameters.  We  also  show  how  to  use  a  ROC  curve  in 
conjunction  with  the  above  relationship  to  determine  an  optimal 
detection  threshold  in  the  signal  processing  system  that  provides 
measurements  to  the  tracker.  A  validation  of  the  modified 
Riccati  equation  is  also  presented. 


1.  INTRODUCTION 


Garden-variety  tracking  problems  involve  processing 
measurements  (e.g.,  range  and  azimuth  observed  by  a  sensor)  from 
a  target  of  interest  and  producing,  at  each  time  step,  an 
estimate  of  the  target's  current  position  and  velocity  vectors. 
Uncertainties  in  the  target  motion  and  in  the  measured  values, 
usually  characterized  as  random  noise,  lead  to  corresponding 
uncertainties  in  the  target  state. 

A  common  and  versatile  approach  to  such  problems  involves 
assuming  that  the  state  dynamics  and  the  measurements  are  both 
corrupted  by  additive,  white,  possibly  Gaussian  noise;  the 
solution  is  then  the  celebrated  Kalman-Bucy  filter  [1-5] ,  which 
is  the  conditional  mean  state  estimator,  best  linear  estimator, 
maximum  &  posteriori  estimator,  maximum  likelihood  estimator,  or 
least-squares  estimator,  depending  upon  one's  point  of  view.  The 
parameters  that  determine  tracking  performance  in  such  a  filter 
are  the  system  matrices  in  the  equations  describing  target  state 
dynamics  and  measurements,  which  will  be  considered  fixed  for  the 
purposes  of  this  discussion,  and  the  covariance  matrices  of  the 
process  and  measurement  noises,  which  specify  the  uncertainties 
in  target  motion  and  measured  values,  respectively. 

In  many  tracking  problems,  particularly  those  arising  in 
surveillance,  there  is  additional  uncertainty  regarding  the 
origin  of  the  received  data,  which  may  (or  may  not)  include 
measurements  from  the  target (s)  of  interest,  interfering  targets, 
or  random  clutter  (false  alarms) .  This  leads  to  the  problem  of 
data  association  or  data  correlation,  which  has  been  attacked  on 
a  number  of  fronts  [6-14]  and  surveyed  in  [15-17].  In  this 
situation,  tracking  performance  depends  not  only  upon  the  noise 
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covariances,  but  also  upon  the  amount  of  uncertainty  in 
measurement  origin.  In  some  of  the  approaches  cited  above 
[6-10] ,  this  dependence  is  explicit  and  is  characterized  in  terms 
of  the  de.tectiQh  BEfltofrility  PD  and  false  alarm  probability  Pp 
(which  is  proportional  to  clutter  density) . 

In  typical  applications,  measurement  data  are  provided  to  a 
tracker  by  upstream  signal  processing  and  detection  algorithms, 
as  indicated  in  Figure  1.  The  process  noise  covariances  are 
normally  selected  on  the  basis  of  experience  and  intuition  (i.e., 
they  are  guessed) .  The  measurement  noise  covariances  are  either 
provided  by  the  signal  processing  algorithm,  as  shown  in  the 
figure,  or  they  are  selected  in  the  same  manner  as  the  process 
noise.  In  any  case,  the  true  noise  levels  are  usually  fixed  by 
target  dynamics  and  sensor  configuration  and  cannot  be  adjusted 
on  line. 

Detection  and  false  alarm  probabilities,  on  the  other  hand, 
are  highly  interdependent  and  adjustable  via  a  detection 
fhisshold:  raising  the  threshold  lowers  both  probabilities,  and 
vice-versa.  This  relationship,  which  also  depends  parametrically 
on  the  signal-to-noise  ratio  (SNR) ,  is  usually  characterized  by 
means  of  a  set  of  receiver  operating  characteristic  (ROC)  curves, 
as  discussed  below  in  Section  4.  The  threshold  is  typically  set 
by  choosing  a  design  point  on  the  applicable  ROC  curve,  based  on 
the  perceived  tradeoffs  between  false  alarms  and  missed 
detections.  However,  to  the  best  of  our  knowledge,  these 
tradeoffs  have  never  included  any  systematic  or  quantitative 
consideration  of  the  effects  downstream  on  data  association  and 
tracking  performance. 

In  this  paper  we  shall  describe  such  a  quantitative 
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FIG.  1.  TRACKING  SYSTEM  BLOCK  DIAGRAM 

relationship.  The  dependence  of  a  tracker's  error  covariance 
upon  the  detection  and  false  alarm  probabilities  is  explicitly 
(but  approximately)  characterized  by  a  scalar  parameter  q2  in  the 
covariance  equation,  called  the  modified  Riccati  equation.  This 
scalar  parameter  depends  upon  the  probabilities  of  detection  and 
false  alarm,  and  also  upon  the  volume  of  the  data  association 
gate,  which  in  turn  depends  on  the  state  error  covariance  matrix 
E.  The  modified  Riccati  equation  can  be  iterated  to  convergence, 
yielding  a  steady-state  £,  and  tracking  performance  can  be 
characterized  by  a  scalar  metric  such  as  determinant  (£)  ,  trace 
(£)  /  or  (in  surveillance  applications)  root-mean-square  position 
error.  This  result  is  important  for  the  following  reasons: 

1.  Contour  plots  of  the  scalar  tracking  performance  metric 
as  a  function  of  detection  probability  and  false  alarm 
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probability  form  a  set  of  tracker  operating 
.character  IP  tic  (TOC)  curves ,  which  can  be  superimposed 
on  ROC  curves  for  the  detector  or  receiver  of  interest 
in  order  to  determine  graphically  the  operating  points 
that  optimize  tracker  performance. 

2.  The  stability  of  the  tracking  process  depends 
critically  on  the  detection  and  false  alarm 

probabilities;  indeed,  a  region  of  apparent  instability 
of  the  modified  Riccati  equation  exists  in  the  PD-PF 
plane  of  the  TOC  curves.  The  implication  of  this  for 
detector/receiver  design  is  that  there  are  settings  of 
the  detection  parameters  that  render  the  output  useless 
for  downstream  tracking. 

3 •  Allocation  of  tracking  resources  (both  computation  and 
communication)  requires  prediction  of  future  state 
error  covariances  under  various  resource 
configurations,  i.e.,  as  a  function  of  detection  and 
false  alarm  probability  and  of  process  and  measurement 
noise  covariance. 

4.  The  same  derivations  provide  a  solution  to  the  related 
problem  of  determining  the  statistical  properties  of 
the  modified  likelihood  function  [18],  used  for 
decision  making  (e.g.  maneuver  detection)  when 
measurement  origins  are  uncertain. 

In  Section  2  the  problem  of  relating  tracking  performance  to 
detection  and  false  alarm  probabilities  is  formulated  in  the 
context  of  probabilistic  data  association.  The  key  TOC  results 
are  derived  in  Section  3,  followed  by  examples  in  Section  4. 
Conclusions  and  suggestions  for  further  research  are  found  in 
Section  5.  Discussion  of  the  numerical  operations, derivation  of 
complex  mathematical  expressions,  and  a  description  of  the 
simulations  used  to  validate  the  modified  Riccati  equation  may  be 
found  in  Appendices  A-E. 
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2.  PROBLEM  FORMULATION 


Consider  a  dynamic  system  (target  model)  of  the  usual  form. 


— k+1  =  £xk  +  0% 


(1) 


yk  =  flxk  +  Zk 


(2) 


where  x  is  the  state  vector,  y  is  the  measurement  vector,  a  and  y 
are  zero-mean,  mutually  independent,  white,  gaussian  noise 
vectors  with  covariance  matrices  Q.  and  £,  respectively,  and  k  is 
a  discrete  time  index.  The  matrices  £,  Q,  H,  Q,  and  £  are 
assumed  known  and  their  dependence  on  k  is  suppressed  here  for 
notational  convenience.  The  initial  state  is  assumed  gaussian 
with  mean  £0 j 0  and  covariance  E0|0.  A  typical  state  vector  would 
include  position  and  velocity  variables,  as  well  as  other 
information  that  relates  to  the  specific  type  of  platform  being 
tracked,  and  typical  dynamics  would  assume  straight-line  or 
great-circle  target  motion  with  disturbances  from  the  process 
noise  &. 

The  tracker's  estimate  of  the  target  state  xk  at  time  k, 
given  data  up  to  time  i,  is  denoted  £kji.  The  error  in  this 
estimate  is  Xk|i  =  Ak-£k|i,  with  error  covariance  matrix  Ek  |  i  = 
E{£k li) ,  where  E  denotes  expectation.  The  discrete-time 
Kalman-Bucy  filter  [2-5]  propagates  these  in  two  stages.  The 
prediction  stage  accounts  for  time  evolution, 


*k|k-l  =  ^k-llk-1 


(3) 
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2k  Ik-!  =  EEk-Hk-lE'  +  SQG' 


(4) 


starting  from  the  initial  conditions  &0|0  and  £0|0.  The  update 
stage  compares  the  incoming  measurement  with  the  predicted 
measurement  £k|k-l  =  H&k|k-1  to  form  the  innovation  vector 

2k  =  ~  |k-l  (5) 

whose  covariance  matrix  is 

ak  4  E{2k2£)  -  HEk  Ifc-ifl.’  +  B  (6) 

The  state  and  covariance  are  then  updated  via 

4  |k  =  *k|k-l  +  (7) 

£k|k  =  +  H^k 

=  Ek|k-1  -  HA**  (8) 

where 

Hk  =  ^klk-lii'Sk1  (9) 

is  the  filter  gain  matrix.  The  resulting  state  estimate,  under 


6 


the  above  assumptions,  is  the  conditional  mean  £k)k  =  E{xk|Yk) 
where  Yk  denotes  all  data  vectors  for  i  i  k. 

Equations  (4)  and  (8)  comprise  the  matrix  Riccatl  equation 
[1-5]  for  the  state  estimation  covariance  that  characterizes  the 
tracker's  performance.  This  iterative  equation  is  deterministic, 
since  for  a  linear  filter,  the  estimation  accuracy  is  independent 
of  the  data. 

In  order  to  avoid  clouding  the  discussion  to  follow,  this 
brief  summary  ignores  a  number  of  complications  that  arise  in 
practice.  If  the  system  is  nonlinear,  for  example,  then  it  can 
usually  be  linearized  and  the  same  basic  equations  can  be  applied 
to  deviations  from  the  nominal  trajectory  [3,  5].  If  the  target 
occasionally  deviates  from  the  assumed  motion  model,  e.g.,  by 
maneuvering,  then  some  decision-making  or  other  machinery  must  be 
provided  to  deal  with  these  instances. 

In  multi-sensor  problems,  the  size  and  composition  of  the 
measurement  vector  often  varies  from  one  time  to  the  next;  in 
other  words,  yk  is  composed  of  independent  subvectors  from 
various  sensors,  any  subset  of  which  may  be  present  at  a  given 
time.  Moreover,  in  the  problem  of  interest  here,  each  sensor 
supplies  not  one  but  several  subvectors  that  must  be  associated 
with  targets.  We  will  avoid  the  resulting  notational  morass  by 
restricting  equations  (5) -(8)  to  apply  to  a  measurement  subvector 
yk  from  a  single  sensor.  In  addition,  j£§.  will  suppress  the  time 
index  k  from  all  variables  except  £,s £,&,  and  Y,  except  where  it 
is  required  for  clarity.  Without  any  loss  of  generality,  the 
data  association  problem  may  now  be  formulated  as  follows. 

At  each  time  step,  the  sensor  provides  a  set  of  candidate 
measurements  to  be  associated  with  targets  (or  rejected).  In 
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most  approaches,  some  preselection  is  done  by  forming  a 
SZfllldatiQn  gate  around  the  predicted  measurement  from  each  target 
and  selecting  those  detections  that  lie  within  the  gate.  There 
are  many  different  approaches  to  establishing  a  correspondence 
between  candidate  measurements  and  targets,  some  of  which  were 
cited  above  in  Section  1. 

In  this  paper  we  shall  focus  on  the  probabilistic  data 
ajsffgcifttibn  (PDA)  method  [6-8,15],  although  the  results  are 
relevant  to  other  methods  [9,  10]  in  which  similar  machinery  is 
used.  The  candidate  measurements  in  a  gate  at  time  fc.  are  denoted 
t  j=l, . . .m,  and  their  corresponding  innovations  are 

2.j  =  -  £/  j=l/ . .  .m  (10) 

The  term  measurement  will  be  used  interchangeably  for  yj  and  £ j , 
since  they  contain  equivalent  information  [5]. 

Considering  a  single  target  independently  of  any  others,  Xj 
denotes  the  event  that  the  j-th  measurement  belongs  to  that 
target  and  Xg  the  event  that  none  of  the  measurements  belongs  to 
it  (no  detection) .  The  PDA  approach  builds  upon  the  assumptions 
that  the  estimation  errors  £  and  £  have  gaussian  densities  at 
each  time  step  (this  is  approximate,  since  there  is  an 
exponentially  growing  tree  of  possible  measurement  sequence 
hypotheses  and  the  true  densities  are  gaussian 
mixtures —  weighted  sums  of  gaussians) .  it  is  also  assumed  that 
the  correct  measurement  is  detected  with  probability  PD 
(independently  at  each  time)  and  that  all  other  measurements  are 
Poisson-distributed  with  parameter  CV,  where  V  is  the  volume  of 
the  validation  gate  and  C  is  the  expected  number  of  false 
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measurements  per  unit  volume.1  with  Note  that  C  =  PF/VC,  where  Vc 
is  the  volume  of  one  resolution  cell  (see  Section  4)  and  Pp  is 
the  probability  of  false  alarm  in  each  cell. 

The  gate  is  normally  a  g-sigma  ellipsoid  g2} 
and  PQ  is  the  probability  that  the  correct  measurement,  if 
detected,  lies  within  the  gate.2  The  gate  volume  is  thus 

V  =  cMgMlak|1/2  (11) 


where  M  is  the  dimension  of  £  and  cM=7rM/2/r(M/2+l)  is  the  volume 
of  the  M-dimensional  unit  sphere  (0-^=2,  c2=ir,  c3=4it/3,  etc.). 

The  conditional  mean  estimate  x  is  obtained  from  (7)  by 
using  the  combined  (weighted)  innovation 

m 

2  s  E  j2j  (12) 

j=l 

where  Pj  =  P{Pj|Yk),  j=0,l,...m,  is  the  posterior  probability 
that  the  j-th  measurement  (or  no  measurement,  for  j=0)  is  the 
correct  one.  These  probabilities  are  given  by  the  following 
expressions  (see  Appendix  C) : 


1Equivalently,  the  number  n  of  false  measurements  has 
probability  mass  function  p(n)  =  e_cv(CV)n/m  and  the  location  of 
each  false  measurement  is  uniformly  distributed  in  the  gate. 

2This  is  just  the  gaussian  probability  mass  in  the  gate,  which 
ls  -P££en  assurosd  to  be  unity  in  practice,  since  Pr.>.99  whenever 
g>M1/2+2. 
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exp(-2«s^12i/2) 

j3j  =  - * - J -  ,  j=l,...m  (13) 

J  m 

b  +  E  exP("5tiak12i/2) 
i=l 

b 

h  =  -  (14) 

m 

b  +  E  expt-SiSk^i/S) 
i=l 

where 

b  =  (27r)M/2C|a|Cl1/2(l-PDPG)/PD 

=  (21r)M/2(CV/cMgM)  (1-PDPQ)/PD  (15) 

The  covariance  update  equation  (8)  is  replaced  by 

Ek Ik  ■  Ek|k-1  -  U-fVHkSuHi;  +  £k  (16) 

where  the  data-dependent  (stochastic)  terms 
m 

Ek  8  Hkl  ZPjSjSJ-JB'JV  (17) 

and  (30  transform  the  original  deterministic  Riccati  equation  into 
a  stochastic  one. 
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dependence  of  tracking  performance,  particularly  the  behavior  of 
^k|k'  on  detecti°n  probability  PD  and  clutter  density  C.  This  is 
accomplished  in  the  next  section  by  means  of  a  deterministic 
approximation  to  the  stochastic  Riccati  equation  (4)  and  (16) . 

Another  major  problem  in  data  association  and  tracking  is 
the  testing  of  hypotheses  for  maneuver  detection,  track 
initiation,  signature  formation,  target  classification,  and  other 
decision-making  purposes.  The  uncertainty  in  measurement  origin 
leads,  in  the  PDA  and  other  approaches,  to  a  modified  likelihood 
function  [18]  involving  the  combined  innovation  2  in  (12) .  A 
major  drawback  of  this  approach  has  been  the  difficulty  of 
computing  the  covariance  matrix  of  2/  but  this  can  now  be  done 
using  an  intermediate  result  to  be  derived  in  the  next  section. 

Finally,  note  that  multiple  targets  can  be  handled 
simultaneously  via  the  joint  probabilistic  data  association 
( JPDA)  approach  [8],  in  which  the  posterior  probabilities 
(13) -(14)  are  computed  jointly  across  a  set  of  potentially 
interfering  targets.  Although  this  is  a  very  important 
extension,  it  greatly  complicates  the  derivations  and  will  not  be 
included  here. 


3.  APPROXIMATE  COVARIANCE  EQUATION 


Since  any  measure  of  tracking  performance  must  depend 
heavily  (or  perhaps  exclusively)  on  the  error  covariance  matrix 
£k|k'  we  sha11  attempt  to  characterize  its  behavior  in  the 
presence  of  uncertainties  in  measurement  origin.  £k|k  is  a 
random  (data-dependent)  matrix  governed  by  the  nonlinear, 
stochastic  difference  equations  (4)  and  (16),  and  hence  its 
behavior  can  only  be  determined  in  a  statistical  sense. 
Moreover,  even  propagation  of  its  first  and  second  moments 
appears  to  be  intractable  except  via  extensive  numerical 
operations. 

Consequently,  we  shall  consider  an  approximation  to  (16)  in 
which  the  random  matrix  £k  defined  in  (17)  and  the  probability  p0 

given  by  (14)  are  replaced  by  their  (prior  to  time  k)  expected 
values 


Ejc  =  E{£k|Yk_1} 

(18) 

P0  =  E{  P0 1  Yk”1}  =  E{P0}  =  1  -  PDPG 

(19) 

where  the  final  expression  is  a  consequence  of  E [P{A|B} ] =P{A} . 

These  substitutions  make  (4)  and  (16)  into  a  set  of 
deterministic  equations  that  can  be  iterated  forward  in  time. 
Because  (16)  is  nonlinear  in  £k|k_lr  this  does  not  yield  E{£k|k}; 
nevertheless,  it  will  give  approximate  values  of  future  state 
error  covariances  in  the  presence  of  uncertain  detections  and 
false  alarms  as  a  function  of  the  environmental  parameters  PD  and 
C,  and  of  the  noise  covariances  £  and  Q. 
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Expansion  of  (18)  yields 


£k  =  E{£k|Yk"1}  =  E{E[£k|m,Yk_1] |yk_1} 

=  £  Et^lmrY^lPCmlY^1}  (20) 

m=0 


where  P{m|Yk_1}  =  P{m}  is  given  by  (44)  in  Appendix  C.  Using 
(17)  and  (12) ,  the  inside  expectation  becomes 


E [£k |m, Yk_1]  = 


Hk  [Hi  (m) -ll2  (m)  ]Hk' ,  m=l,2. 


0, 


m=0 


(21) 


where 


m 

U-idn)  =  E[  £  0^21 1  m^”1] 
i=l 

and^ 


m  m 

a2(m)  =  E[£  Pi2iE  PjS^lmrY^1] 

i=l  j=l 


(22) 


3 The  second  expression  in  (23)  is  obtained  as  a  intermediate 
result  in  Appendix  D. 
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(23) 


m 


=  E[E  Pi2i2ilmfYk_;L 

i=l 


1 


The  expected  values  are  obtained  by  multiplying  the  quantities  in 
square  brackets  by  the  joint  prior  density  p(£j , . . .2m|m, Yk“^) 
given  in  (49)  and  integrating  over  the  validation  gate. 

Considerable  simplifications  result  if  one  applies  a  linear 
transformation  of  variables  (S^^2  )to  obtain  a  spherical  gate, 
followed  by  a  change  to  spherical  coordinates.  Because  of  the 
spherical  symmetry  of  the  gaussian  density  and  of  the  expressions 
(22) -(23),  off-diagonal  elements,  cross-terms,  and  angular 
variables  drop  out  like  flies,  leaving  scalar  integrals  over  the 
radial  variables  (i.e.,  over  ||£j||,  j=l,...m).  The  detailed 
derivation,  which  is  given  in  Appendix  D,  leads  to 


U.!  (m) 
U2(m) 


m 


'M 


PDPGm  +  (1-PDPQ)CV  (2tt)m/2 


^k 


m 


'M 


PDPQm  +  (1-PDPG)CV  (21 r) 


M/2 


(-. )“ 


I2(m)-k 


(24) 


(25) 


where  is  the  covariance  matrix  of  the  correct  innovation  and 
the  scalar  integrals  1^  and  I2(in)  are  defined  as 


=  H  rM+1exp(-r2/2)dr 


(26) 
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g  g  expc-rfjr!2 

*2  (m)  f  •••  r  "  ~ 1  ■  ■  try"  —  \  M— 1  -j 

\  l  _  ^rlr2***rm^  dr^...dr 


m 


0  b  +  2Z  exp(-r-j2/2) 

j=l 


m 
(27) 


and  the  constant  b  =  (2tt)  M/2C  |ak  | 1/2  (1-PDPQ) /PD  was  defined  in 
(15)  . 

Substituting  (24) -(27)  and  (44)  into  (20)  and  cancelling 
leads  to 


-  (q1-q2)HkakHl; 


(28) 


where 


.  „  Jqj _  «  ^icv)1"-1  _  cM 

^  D  (2tt)m/2  m=l  (m-1)  I  Zl  Pd  (2tt)M/2  Ix 


(29) 


CM  00  e“CV(CV)m“1  /  M\m-1 

- 7TF>  £  - l  —  )  i2(®> 

(2tt)m/2  m=l  (m-1)  1  \gM/ 


(30) 


and  substitution  of  this  and  (19)  into  (16)  yields  the 
deterministic  equation 


^k|k  £k|k-l  "  ^PDPG“ql+q2^k^kHk  (31) 

This  can  be  simplified  further  by  noting  that  for  typical 
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values  of  g  and  M  (g=4  or  5  and  M<10)  ,  Pq  is  approximately  1  and 
<3l  is  approximately  PD  (substitute  x1/2  for  r  in  ij.  to  get  a 
gamma  function) . 

The  upshot  of  all  this  is  that  the  deterministic 
approximation  to  the  covariance  equation  becomes 


^k|k  ^k|k-l  “  ^^k^k^lc  (32) 

where  the  scalar  q2  lies  between  0  and  1  and  depends  on  P^,  c, 
and  V,  the  volume  of  the  validation  region  at  time  k.  Comparing 
this  to  (8) ,  it  is  clear  that  the  factor  q2  reduces  the 
covariance  improvement  due  to  the  term  wsw 1 ;  the  smaller  q2  is, 
the  greater  the  degradation. 

Since  C  is  proportional  to  Pp  and  V=cMgM |£k | x/2  from  (11), 
equations  (4)  and  (32)  may  be  written  as 

Ek|k-1  =  Ok-Hk-iE'  +  fiQS' 

(33) 

^k|k  ^k|k-l  "  ^2  (S-k;PDf  pF^kSk^lc 


where  Hk  and  £k  depend  upon  £kjk-1  via  (6)  and  (9).  This 
modified  Riccntj  equation  describes  (approximately)  the  behavior 
of  the  PDA  tracking  filter  as  a  function  of  the  detection  and 
false  alarm  probabilities  P^  and  Pp.  The  approximation  is 
validated  via  Monte  Carlo  simulations  in  Appendix  E.  We  shall  now 
use  it  to  characterize  the  dependence  of  tracking  performance  on 
PD  and  Pp. 
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For  most  values  of  PD  and  Pp,  (33)  can  be  iterated  until  it 
converges  to  a  steady-state  covariance  matrix  £(PD,Pp)  (the 
stability  issue  is  discussed  below) .  In  order  to  obtain  a  scalar 
tracking  performance  metric,  one  can  then  extract  the  steady- 
state  root-mean-square  (RMS)  position  error 


(34) 


where  p^  and  p 22  are  the  diagonal  elements  of  £  that  correspond 
to  target  position. 

We  shall  refer  to  a  contour  plot  of  (34)  as  a  tracker 
operating  characteristic  (TOC) .  This  name  is  chosen  because  the 
well-known  Lsgeiver  operating  characteristic  (ROC)  curve  in  the 
same  P[)“Pp  Plane  is  the  locus  of  possible  operating  points  for  a 
detector/receiver,  where  a  particular  operating  point  on  the 
curve  is  determined  by  the  detection  threshold  level.  Thus,  if 
the  ROC  curve  is  superimposed  on  the  TOC  contours,  the  dependence 
of  tracking  performance  on  detection  threshold  can  be  determined 
directly.  This  will  be  illustrated  in  the  next  section  with  an 
example. 

There  are  various  other  performance  metrics  that  can  be 
used,  of  course,  such  as  the  determinant  or  trace  of  £.  In  many 
applications,  the  steady-state  covariance  may  not  be  appropriate: 
one  can  instead  use  the  value  of  £k|k  obtained  by  iterating  (33) 
over  a  fixed  period  of  time  from  a  standard  £0  j  0 . 
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4.  EXAMPLES 


The  target/sensor  geometry  shown  in  Figure  2  was  used  in  the 
multi-target  tracking  examples  of  [8].  Taking  the  (linearized) 
values  of  £,  £,  fl,  Q,  and  £  from  the  initial  time  in  that 
example,  we  have  iterated  (33)  to  obtain  the  steady-state  RMS 
position  error  (34)  for  various  values  of  PD  and  Pp.  Evaluation 
of  <32 (Sk'pD' pf)  was  carried  out  using  a  look-up  procedure  from 
tables  generated  off-line  (see  Appendix  A  for  details) . 

Tracker  operating  characteristics  will  be  shown  for  two 
different  measurement  types,  in  the  first  example,  the  target  is 
tracked  using  measurements  of  bearing  (azimuth)  and  frequency 
from  sensors  20  and  22  at  5-minute  intervals.  The  process  noise 
matrix  (GOG 1 )  is  diagonal,  with  standard  deviations  of  .2°  in 
course,  .2  knots  in  speed,  and  .01  Hz  in  source  frequency. 

The  measurement  noise  matrix,  also  diagonal,  has  a  standard 
deviation  of  5°  in  bearing  and  .08  Hz  in  frequency.  We  further 
assume  that  the  sensor  signal  processing  is  able  to  resolve 
signals  separated  by  about  4°  and  .15  Hz.  We  may  thus  view  the 
space  of  bearing/frequency  measurements  as  a  collection  of 
resolution  cells,  with  the  tracker's  validation  gate  encompassing 
some  subset  of  these  cells.  in  practice,  the  detector/receiver 
will  have  an  ad  hoc  rule  prohibiting  detections  in  adjacent 
cells,  so  that  the  effective  cell  volume  is  about  Vc  =  ,3°-Hz. 
We  assume  further  that  false  alarms  occur  independently  in  each 
cell  with  probability  Pp,  so  that  the  clutter  density  (expected 
number  of  false  alarms  per  unit  volume)  is 


C  =  Pp/.3°-Hz 


(35) 


With  these  assumptions,  the  TOC  contours  shown  in  Figure 
3  have  been  computed  (note  that  Pp  ranges  only  from  0  to  0.1). 
This  is  a  contour  plot  of  (34),  with  performance  improving  (i.e., 
position  error  decreasing)  toward  the  upper  left-hand  corner. 
Performance  degrades  in  the  other  direction,  as  PD  decreases 
and/or  Pp  increases,  and  there  is  a  region  in  which  the  modified 
Riccati  equation  (33)  does  not  appear  to  converge  to  a  finite 
steady-state  covariance  £.  The  question  of  stability  is 
discussed  further  in  Appendix  B. 

Figure  3  specifies  tracking  performance  as  a  function  of  PD 
and  Pp.  in  order  to  determine  what  values  of  these  probabilities 
are  achievable,  we  need  LScejyqr  operating  characteristic  (ROC) 
curves  for  the  detection  system  that  provides  measurements  to  the 
tracker.  To  this  end,  we  shall  assume  that  the  detection 
algorithm  is  equivalent  to  a  set  of  classic  quadrature  receivers 
or  incoherent  patched  filters  [19] ,  one  operating  on  each 
resolution  cell  in  bearing/frequency  space.  The  quadrature 
receiver  assumes  a  sinusoidal  signal  of  unknown  phase.  Under  the 
signal-plus-noise  hypothesis,  the  test  statistic  has  a  Rician 
distribution,  which  reduces  to  a  Rayleigh  distribution  in  the 
noise-only  case.  Expressions  for  PD  and  Pp  may  be  derived  [19] 
and  used  to  compute  the  ROC  curves  shown  in  Figure  4. 

For  a  given  signal-to-noise  ratio  (SNR) ,  the  corresponding 
ROC  curve  is  the  locus  of  possible  operating  points  that  the 
detector  can  assume,  depending  on  where  one  sets  the  detection 
threshold.  In  Figure  3  one  such  curve  (SNR=8  dB)  is  superimposed 
as  a  dashed  line  on  the  TOC  contours.  This  shows  graphically  how 
tracking  performance  depends  on  the  operating  point  of  the 
detector  (i.e.,  on  the  detection  threshold).  in  particular, 
performance  is  optimal  at  the  operating  point  indicated  by  a. 
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There  is  a  relatively  broad  region  about  this  point  where 
performance  is  near-optimal,  but  performance  degrades 
significantly  thereafter. 

Coherence  measurements 

In  the  second  example,  the  target  is  tracked  by  cross- 
correlating  signals  between  pairs  of  sensors  to  obtain 
measurements  of  time  delay  difference  and  Doppler  difference  from 
sensor  pairs  20/21  and  21/22  (see  Figure  2)  at  5-minute 
intervals.  Standard  deviations  of  4  sec  in  time  difference  and 
.004  Hz  in  Doppler  difference  are  assumed,  and  the  effective 
resolution  cell  volume  is  .008  sec-Hz,  so  that 


C  =  Pp/.008  sec-Hz 


(36) 


This  leads  to  the  TOC  contours  shown  in  Figure  5. 

Analysis  of  the  cross-correlation  algorithm  is  somewhat  more 
difficult  than  that  of  the  quadrature  receiver.  Nevertheless,  if 
we  assume  that  both  signal  and  noise  are  sample  functions  from 
white,  gaussian,  random  processes  and  that  the  time-bandwidth 
product  is  500  sec  x  .25  Hz,  we  can  obtain  the  ROC  curves  shown 
in  Figure  6.  The  parameter  is  now  coherence  between  the  two 
channels,  rather  than  SNR. 

Again,  for  a  given  coherence,  the  corresponding  ROC  curve  is 
the  locus  of  possible  operating  points  for  the  coherence 
detector.  In  Figure  5,  the  curve  corresponding  to  a  coherence  of 
.025  is  superimposed  as  a  dashed  line  on  the  TOC  contours  to  show 
graphically  how  tracking  performance  depends  on  the  operating 
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point  of  the  coherence  detector  (i.e.,  on  the  detection 
threshold) ,  and  the  optimal  point  is  indicated  by  a. 
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FIG.  2.  TARGET-SENSOR  GEOMETRY 
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FIG.  3.  TRACKER  OPERATING  CHARACTERISTIC  (TOC)  CONTOURS  FOR 

BEARING/ FREQUENCY  MEASUREMENTS  (DASHED  LINE  IS  SUPERIMPOSED 
ROC  CURVE  FROM  FIGURE  5) 


FIG.  4.  ROC  CURVES  FOR  QUADRATURE  RECEIVER 
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TRACKER  OPERATING  CHARACTERISTIC  (TOC)  CONTOURS  FOR 
TIME/DOPPLER  DIFFERENCE  MEASUREMENTS  (DASHED  LINE  IS  SUPERIMPOSE 


FIG.  6.  ROC  CURVES  FOR  COHERENCE  RECEIVER 
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5.  CONCLUSION 


We  have  established,  for  the  first  time,  an  important 
relationship  between  thresholds  in  detector/receivers  and 
performance  in  downstream  trackers.  More  specifically,  a 
modified  Riccati  equation  determines  the  approximate  state  error 
covariance  of  a  probabilistic  data  association  (PDA)  tracking 
filter  as  a  function  of  the  threshold-dependent  probabilities  of 
detection  and  false  alarm.  By  plotting  contours  of  tracking 
performance  (in  this  case,  steady-state  RMS  position  error)  in 
the  P£)“Pp  plane  and  then  superimposing  a  ROC  curve  for  a 
particular  SNR,  one  can  determine  the  optimal  detection  threshold 
graphically. 

Several  extensions  of  this  concept  are  of  interest.  The 
graphical  method  for  selecting  an  operating  point  can  be  replaced 
by  a  mathematical  optimization:  an  obvious  necessary  condition  is 
that  the  ROC  and  TOC  curves  be  tangent.  However,  the  practical 
difficulty  of  computing  the  required  differentials  to  solve  the 
necessary  conditions  is  substantial.  An  approximate  (e.g.  table 
look-up)  procedure  for  optimization  would  be  useful  for  dealing 
with  multi-dimensional  TOCs,  such  as  will  occur  if  different 
receivers  are  allowed  to  have  different  thresholds  or  if  bearing/ 
frequency  and  time/Doppler  measurements  are  used  simultaneously. 

Another  issue  of  major  importance  is  the  optimization  of 
tracking  performance  when  the  signal's  SNR  is  not  known.  In  this 
case,  several  ROC  curves  are  involved  and  one  must  select  a 
threshold  that  is  best  (in  some  sense)  for  a  whole  range  of  SNRs. 
Alternatively,  an  adaptive  thresholding  scheme  can  be  devised, 
whereby  the  SNR  is  monitored  and  the  threshold  adjusted  so  as  to 
maximize  performance  along  the  current  ROC  curve. 
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The  results  obtained  here  apply  to  the  probabilistic  data 
association  framework,  in  which  each  target  is  considered 
individually  in  the  presence  of  random  clutter  (false  alarms) . 
It  will  be  useful  to  modify  these  results  to  deal  with  other  data 
association  schemes,  particularly  the  joint  PDA  approach  [8],  in 
which  multiple  interfering  targets  are  accounted  for  by  computing 
the  posterior  probabilities  (13) -(14)  jointly  across  a  set  of 
targets. 

A  Monte— Carlo  test  of  the  validity  of  the  key  approximation, 
in  which  £k  and  were  replaced  in  (16)  by  their  expected  values 
is  presented  in  Appendix  E,  thus  validating  the  entire  ROC/TOC 
approach. 
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APPENDIX  A 

PROPAGATION  OP  THE  MODIFIED  RICCATI  EQUATION 


Propagation  of  the  modified  Riccati  equation  (33)  requires 
evaluation  of  the  scalar  parameter  q2r  defined  by  (27)  and  (30) . 
This  was  carried  out  numerically  as  follows.  First,  a  50,000- 
point  Monte-Carlo  integration  scheme  was  used  repeatedly  to 
create  a  table  of  values  of  the  integral  (27)  for  various  values 
of  b,  m=l-15,  and  M=l-4  (the  gate  size  was  fixed  arbitrarily  at 
g=4).  Then  routines  from  the  numerical  package  IMSL  were  used  to 
compute  spline  coefficients  for  interpolation  over  b.  Finally, 
the  tables  and  coefficients  were  embedded  in  a  subroutine  that 
accepts  values  of  M,  PD  and  CV  (=PFV/VC)  and  evaluates  (30)  by 
truncating  the  sum  at  m=15  and  using  the  spline  interpolation  to 
evaluate  the  integral. 

As  can  be  seen  from  (27)  and  (30) ,  for  a  given  dimension  of 
the  measurement  vector,  M,  and  gate  size,  g,  the  resulting 
coefficient  q2  depends  only  on  the  target  detection  probability, 
PD,  and  the  expected  number  of  false  measurements  in  the  gate, 
CV.  Thus  the  results  obtained  from  the  above  evaluation  yield 
"universal  curves",  presented  in  Figure  7. 

With  this,  propagation  of  (33)  was  straightforward. 
However,  because  this  form  of  the  right-hand  side  of  the  update 
equation  is  known  to  have  poor  numerical  properties,  it  was 
replaced  in  the  computer  program  by  the  equivalent  equation 


*k|k  -  d-^klk-l  +  q2l<i-“kfi>2k|k-l<i-HkH>’  +  (37) 
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To  prove  the  equivalence,  one  may  simply  use  the  well-known 
identity  given  by  (8)  on  the  terms  in  square  brackets. 


CV 

Figure  7.  Factor  q2  for  measurement  of  dimension  M=2  and  gate  g=4 
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APPENDIX  B 

STABILITY  OP  THE  MODIFIED  RICCATI  EQUATION 


The  standard  Riccati  equation  defined  by  (4)  and  (8)  is 
well-known  to  be  stable,  provided  appropriate  controllability  and 
observability  conditions  hold  [3];  it  will  always  converge  to  a 
steady-state  covariance  £.  The  modified  Riccati  equation  (33) , 
with  an  additional  covariance-dependent  factor  q2  on  the  right- 
hand  side,  presents  some  very  challenging  and  unresolved 
stability  issues.  It  is  clear  from  our  numerical  results  that 
(33)  converges  to  a  steady-state  £(PD,Pp)  for  most  values  of  the 
parameters  PD  and  Pp.  In  the  instability  regions  indicated  in 
Figures  3  and  5,  (33)  diverges  numerically,^  but  this  does  not 
necessarily  imply  that  the  equation  is  unstable  in  a  mathematical 
sense. 


Although  we  have  as  yet  been  unable  to  establish  any  general 
theorems  on  either  stability  or  instability,  examination  of  a 
scalar  example  provides  some  useful  insights. 

Suppose  that  x  and  y  are  both  1-dimensional  with  F=G=H=1  and 
Q,R>0,  in  which  case  (33)  reduces  to 

pk+l|k  =  pk|k-l  "  S2pk |k-l/tpk |k-l+Rl  +  Q  (38> 

Note  that  F=1  corresponds  to  an  integrator,  so  that  the 
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Once  q2  underflows,  divergence  is  inevitable. 
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covariance  diverges  unless  it  is  updated.5  if  (38)  converges  to  a 
steady-state  value  P,  then  we  can  substitute  Pk+1|k  =  pk  |  k-1  =  p 
to  obtain  the  quadratic  equation 

q2P2  -  QP  -  QR  =  0  (39) 


which  has  a  positive  solution 


P(q2) 


Q  +  Q2+4QRq2 
2q2 


(40) 


It  can  be  shown  that  for  any  fixed  value  of  q2  between  0  and  1, 
(38)  is  stable  and  converges  to  (40)  from  any  initial  Pn0>0. 
Indeed,  we  conjecture  that  (33)  is  stable  for  any  fixed  value  of 
q2/  and  that  instabilities  occur  only  because  of  the  dependence 
of  q2  on  P. 

The  stability/instability  question  can  be  illuminated  for 
this  scalar  example  by  plotting  P  vs.  q2.  in  Figure  8,  the  curve 
labeled  P(q2)  has  been  computed  from  (40),  while  the  curve 
labeled  q2(P)  has  been  evaluated  from  (30)  using  the  subroutine 
described  above  in  Appendix  A,  for  PD=.7  and  Pp=.l.  One  approach 
to  finding  a  steady-state  solution  (if  one  exists)  of  (38)  is  to 
pick  an  initial  P>0  and  then  alternately  evaluate  q2(P)  and 
P(q2)f  as  indicated  by  the  arrowed  paths.  This  procedure  is 


5 

In  practically  all  tracking  problems,  the  plant  equation 
contains  at  least  one  integration  from  velocity  to  position. 
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stable  for  P<800,  in  which  case  it  converges  to  the  lower 
intersection  of  the  curves.  For  P>800,  the  procedure  diverges  as 
indicated,  suggesting  that  (38)  is  unstable  for  large  enough 

P1 1 0  * 

Figure  9  shows  the  same  P(q2)  curve  and  a  new  q2(P)  curve, 
computed  for  Pjj=.2  and  Pp=.l.  In  this  case  there  is  no  point  of 
intersection  and  the  arrowed  pathways  diverge  for  any  initial 
value  of  P. 

These  examples  and  others  we  have  seen  strongly  suggest  that 
a  stable  intersection  is  always  accompanied  by  an  unstable  one  at 
a  higher  value  of  P  if  Pp>0.  This  implies  that  even  a  small 
amount  of  clutter  will  render  the  equation  unstable  for 
sufficiently  large  P.  We  conjecture  that  the  same  statement  holds 
for  the  multivariable  case  of  (33) . 
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FIG.  8.  STABLE  AND  UNSTABLE  INTERSECTIONS 


FIG.  9.  UNSTABLE  SITUATION 
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APPENDIX  C 

PROBABILITY  CALCULATIONS 


In  this  appendix  we  shall  derive  a  number  of  expressions 
needed  in  the  main  text.  Letting  7j (m)  denote  the  prior 
probability  of  the  event  Xj*  conditioned  on  m,  the  total 
probability  theorem  yields 


Yj(m)  =  P{Xj|mrYk_1}  =  P{X j |m} 


=  P{x j |mF=m-l,m}P{mF=m-l |m}  +P{ Xj |mF=m,m}P{mF=m |m} 


■{ 


(1/m) P{mF=m-l |m}  +  (0) P{mF=m|m} , 
(0) P{mF=ra-l |m}  +  (1) P{mF=m |m} , 


j “1 , • .  •  m 
j=0' 


(41) 


because  mF,  the  number  of  false  measurements,  must  be  either  m-1 
(if  the  target  is  detected)  or  m  (if  it  is  not).  Using  Bayes' 
rule  and  the  assumed  Poisson  distribution  for  false  measurements, 


P{mF=m-l |m} 


P { m | mF=m-l } P { mF=m-l } /P { m } 

[PDPG] [e-cv (CV) m~l/ (m-1) l]/P{m) 
PpPGm/ [  PjjPGm  +  (1-PDPG)CV] 


(42) 


32 


P{mF=m|m}  =  P{m |mF=m}P{mF=m}/P{m} 

=  [1-PDPG] [e_cv(CV)m/m!]/P{ra} 

“  (l*-PDPG)CV/tPDPGm  Pj)PG)CV]  (43) 


where  the  denominator  P{m}  is  the  prior  probability  of  m  and  is 
equal  to  the  sum  of  the  numerators  in  the  two  equations: 

P{m}  =  P{m|Yk_1} 

=  [PDPGm  +  (1-PDPG) CV] e_CV (CV) m-1/m! ,  m=0,l,...  (44) 


Substituting  back  into  (41)  yields 


pDpG/[PDpGm  +  <1-PDPG)CV], 
(l-PDPG)CV/[PDPGm  + 


j=l, . . .m 
j=0 


(45) 


Note  that  Yj (m)  is  independent  of  j  for  j>0. 

Using  Bayes'  rule,  the  posterior  probabilities  in  (12)  can 
be  expressed  as 


Pj  P{Xj|Y  }  —  P{  Xj  12^  /  •  •  ^ } 

"  P(2l»- --2m lx j.mfYk"1)p{Xj |m, Yk_1}/p(2lf . . .2m|mrYk_1) 
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(46) 


=  P(£ir • • .2ml Xj fmfyk-1) 7j (m)/  £  numerators 

j=0 


The  first  factor  is  the  joint  probability  density  of  the  m 
candidate  measurements ,  conditioned  on  the  j-th  one  being 
correct.  According  to  the  PDA  assumptions,  the  correct 

measurement  has  a  gaussian  density 

N(£j;fi,Sk)/PG  =  (l/PG)exp(-£<S_12y2)/(2Tr)M/2|£kl1/2  (47) 

with  mean  a  and  covariance  S^,  where  the  factor  1/Pq  accounts  for 
its  restriction  to  the  validation  gate,  and  each  incorrect 
measurement  has  a  uniform  density  V”1.  It  follows  that 


P(2l, 


.£mIXj/m,Yk"1) 


v'm+1N(2j;a,ak)/pG, 

>m, 


j=lf 

j=0 


m 

(48) 


The  second  factor  in  (46)  is  the  prior  probability  of  Xj,  given 
by  (45)  .  The  denominator  is  the  joint  prior  density  of  the 
measurements,  conditioned  only  on  m  (and  the  past  data) , 


P(21'--*S:nilm,yk~;L)  * 


V“%(m)  +  V"m+1  2  (l/PG)N(2j;£,&k)Y.i(iii) 

3=1 


m 


(49) 
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Note  that  with  the  above  conditioning,  the  validated  measurements 
are  not  independent,  i.e.,  (49)  is  not  equal  to  the  product  over 
j  of  the  marginal  Rtlpr,  densities 


P(2j  |m,Yk-1)  =  V’^l-Kj  (m)  ]  +  (1/PG) N (2j  7j  (m)  (50) 

Finally,  substitution  of  (45),  (48),  and  (49)  into  (46) 

followed  by  a  certain  amount  of  rearrangement  yields 


exp(-2»sJ12j/2) 

S-j  ”  — ^  j— 

J  m 

b  +  £  exp(-2«S^12i/2) 
i=l 


(51) 


00 


b 


b  + 


m 


E 

i=l 


exp(-£«ak12i/2) 


where 

b  =  (2ir)M/2C|akl1/2(l-PDPG)/PD 

=  (2tr)M/2(CV/cMgM)  (1-PDPG)/PD 


(52) 


(53) 
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APPENDIX  D- 

DERIVATION  OP  THE  U±  AND  U2  INTEGRALS 


In  this  appendix  the  expectations 

yi(m)  £  e[  |  ei2i2ilm>Yk’1] 


(54) 


and 


U2(m) 


(55) 


will  be  evaluated,  where  is  given  by  (13)- (14)  and  the  joint 
density  of  y^/...ym  is  (48).  In  order  to  simplify  the  arguments, 
we  shall  make  use  of  the  fact  that  S  is  positive  definite  and 

J\ 

change  variables  to 


(56) 


so  that 


t  r*  — 1~  _  /  w 

-i-k  -i  -  EiEi 


(57) 


and 


y  .y  .r,'. 

iiii  -k-i-inc 


(58) 
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In  terms  of  the  new  variables  the  validation  gate  becomes 
a  sphere  U  :  ||.£  ||  <  g  1  with  volume  cMg  =  V;  the  can 
be  rewritten  as 


A 


r*  ll  Sill 


m  777:  it  2  ' 


b  +  l  e”*5  H  -j  II 
j-1 


i  =  1, . . .m 


(59) 


where 


b^  (2iT)M/2(CV/cMgM)  (1-PDPG)/PD 


(60) 


is  the  same  as  in  (15)  .  Note  that  the  dimensionless  quantity  CV, 
the  average  number  of  false  measurements  in  the  (now  spherical) 
gate,  is  unaffected  by  the  variable  change.  The  joint  density 
(49)  becomes 


P  ( ^  2. ' 


=  Y  o (m)V_m  +  Yj(m) V~m+1 


m 

l  mi*  ;0,I)/Pr(61) 

i=l  -1  -  -  '7 


if  ||  5±||  £  g  for  all  i,  and  0  otherwise. 
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(21) 


Using  the  change  of  coordinates  from  y  to  expressions 
through  (23)  may  be  reexpressed  as  follows: 


k-1. 


E[Pjm,Y"  XI  = 


U1(m) 


02  (m) 


m 


E[i=i§i-i  = 


(62) 


(63) 


(64) 


Let  D  stand  for  the  spherical  validation  gate  region.  Then  the 
above  two  expectations  may  be  written  out  explicitly  as  multiple 
integrals  of  certain  matrices,  taken  over  m  copies  of  D: 


Hi  - 


....  J  A(£^, . . .  r£m)  p  ,  *»m)  <“^1*  *  »dCj 

D  D  d 


m 


(65) 


H2  = 


D  D 


— lf  *  *  * ' -m^  P ^ i ,  •  •  •  ,  ^Hi *  *  * ( 66 ) 


where  p  is  the  joint  innovations  density  (61)  with  the  conditioning 


pq 


B 


and  A,B  are  Mxm 

matrices  described 

componentwise  by 

m 

j,  t&i 

1  £  p,q  <_  M 

(67) 

m 

1,3=1  1  3  1  3 

1  £  p,q  £  M 

(68) 
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Here  the  M-dimensional  vector  has  been  denoted  as  (S J, ,t 
The  integrals  (65)  and  (66)  are  actually  mM-fold  integrals. 
Although  the  complexity  of  undertaking  such  a  large  computation 
directly  would  prove  formidable,  a  number  of  observations  show 
that  these  integrals  have  a  fairly  simple  and  straightforward 
structure: 

Observati on  1 .  The  matrices  6^  and  U2  are  diagonal. 

From  the  expressions  (67)  and  (68)  one  may  deduce  that 
off-diagonal  elements  will  integrate  out  to  zero.  For  when 
P  7*  q,  both  A  and  B  become  odd  polynomials  of  second  degree 
in  the  variables  ^,...,5^,  with  coefficients  either  single 
8's  or  pairwise  products  of  0's.  Now  both  the  §'s  and  the 
two  terms  of  the  probability  density  p  in  (61)  are  positive 
functions  depending  only  on  \\q^\  , . . .  ,||  Sm||  ,  i.e.  they  have 
the  same  values  at  antipodal  points  of  a  sphere.  The  odd  - 
polynomials  from  A  and  B  will  have  opposite  values  at 
antipodal  points,  so  that  their  contributions  to  the  total 
integral  will  cancel. 

Physically,  this  amounts  to  the  observation  that  off  -diagonal 
elements  of  the  inertia  tensor  vanish  in  a  principal  axis  system, 
for  a  spherical  body  with  the  shape  of  the  validation  region  and 
with  mass  distribution  density  appropriately  defined  in  terms  of 
the  0's  and  p. 
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Observation  2 .  Integrals  of  cross  terms  (i^j)  vanish  in  the 
sum  (68)  defining  B 

y  pq 

From  the  previous  observation,  we  need  only  concern 
ourselves  with  diagonal  entries  B  .  Again  when  i?*j  we  have  an 
odd  polynomial  in  the  C's,  and  the  same  arguments  as  before 
apply  again  to  show  its  integral  over  a  sphere,  weighted  by 
the  density  p,  must  vanish. 

From  the  point  of  view  of  the  original  definition  (64) 

A 

of  Uj*  we  could  say  that  distinct  innovations  are  orthogonal: 

their  inner  products  vanish  and  do  not  contribute  to  the  second 
moment  of  the  combined  innovations. 


Observation  3.  Each  term  of  A  or  B  making  a  nonvanishing 
contribution  to  U^,  N=1  or  2,  has  the  same  value,  given  explicitly 
by 


a 

» 

e-1'  1 

j 

0 

b+l  e~hri2_ 

3-1 


(69) 


where 


C  =  C  (m,M)  = 


(2tt  ) 


Yj  (m) 

M/2  ,_  _Mxm-l„ 


(cM9  I 


(70) 
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(71) 


A  Y0 (m) 

C  =  C  (m,M)  &  — °-  „  m 

u  U  icH<P)m 


K  -  K(m,M)  =  Mm"1c2 

M 


(72) 


To  prove  this,  let  us  consider  the  first  term  of  the  (1,1) 
entry  of  each  integral,  which  can  be  written  as  6^U*) 2p , . . . , 5  ) , 
N=l,2.  Since  p  and  depend  only  on  the  norms  ||?i||  =  r.^, 

it  is  clearly  advantageous  to  introduce  spherical  coordinates 
in  each  of  the  m  copies  of  the  5— space: 


1 2 


M 


d£i  =  d;tdcf...d;?  =  r^"1dridm“ 


M-l 


(73) 


M-l 


Here  da^  denotes  the  "unit  solid  angle,"  i.e.  the  surface 
element  on  the  unit  sphere  (S^)2  +  (£2)2  +  ...  +  (51?)  2  =  1. 


Making  these  substitutions,  the  integral  for  the  expected 


value  of  8^(c;^)2  becomes 


a 

f 


e 


m  .  2 
|b+  l 

j=l  J 


(?i)2<cgjI1e”5r3+cu)  <rridrid“51 


m  -%r2 


M-l 


M-l.  ,  M-l ,  _  M-l 

)  . . .  ( r  dr  dco 

m  mm 


(74) 
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The  constants  C^C^,  and  K  defined  in  (70)  ,  (71)  are 

introduced  to  simplify  the  notation. 


The  expression  in  (74)  is  not  yet  complete  because  the 
variable  still  has  to  be  reexpressed  in  angular  co^-coordinates 
for  the  first  copy  of  the  M-dimensional  state  space.  However, 
the  other  angular  variables  ok?”1,  . . .  ,a)M~1  do  not  occur  in  the 

^  III 

integrand,  so  that  their  integrals  can  be  absorbed  in  the 
constant  K.  The  fact  needed  to  do  this  is  that  the  integral 

of  the  unit  solid  angle  is  just  the- surface  area  of  the  unit 
sphere,  which  is  M  times  c^  (the  volume  of  the  unit  sphere) . 


Figure  10.  Unit  sphere  SM 


To  handle  the  £  ^-dependence 
of  the  integrand  in  (74)  we  are 
going  to  use  an  argument  almost 
identical  to  the  one  used  to  derive 
the  volume  of  the  unit  sphere 
(which  is  to  start  from  a  known 
integral  and  "work  backwards") . 

The  known  integral  that  we  wish  to 
exploit  is 


1  9  M 

(?7)  exp  (-%  l 

j-1 


2  12 
Cj)dCJdt£. 


.dC“  =  (2tt)M/2 


(75) 
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Now  if  this  integral  is  rewritten  in  spherical  coordinates,  we 
should  make  a  substitution  C*  =  r^  (co^"1)  where  f  is  a  smooth 
function  independent  of  r^.  For  example  if  coordinates  are 
chosen  on  the  sphere  as  depicted  in  Figure  10 ,  we  would  have 

)  =  cos  0,  where  the  angle  9  is  measured  downward  from  the 
north  pole  as  shown.  Rewriting  (75)  in  spherical  coordinates 
gives 


<2*>M/2  =  }  rf1  dr, 

0 


f2«of1)d<»»-1 

sphere 

rlEl 


(76) 


(75) 

gamma 


2 

The  substitution  Hx1  =  x1  transforms  the  radial  part  of 
into  the  well-known  integral  representation  for  the 
function 


r  (t) 


(77) 


so  now  we  have 


>M/2  0M/2„  ,M  .  .  .  f  c2  .  M-l.  ,  M-l 

(2tt)  =  2  '  T  (^  +  1)  f  (<u1  x)dco1  (78) 

sphere 

rl=1 

We  are  now  in  a  position  to  solve  for  the  nonradial  part  of  the 
integral  in  (78)  .  But  when  this  is  done,  the  quantity  we 
obtain  is  recognized  to  be  cM,  the  volume  of  the  unit  sphere: 
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Actually,  by  using  properties  of  the  gamma  function,  this 
expression  for  c^  can  be  written  more  explicitly  as 


2M/27rM/2 

M(M-2)  .  .  .4*2 


M/2 

IT 

if  M  is  even 

(M/2)  ! 


2tt 

M  CM-2 

2  (M+l) /2  ^(M-D/2 

M (M-2) .. .5‘3*1 

l 


if  M  is  odd 


(80) 


Now  returning  to  the  original  integral  (74)  ,  we  note  that 
even  though  the  integrand  is  not  the  same  as  (75)  ,  the  angular 

dependence  is ,  so  that  (79)  may  be  applied.  Absorbing  this 
cfl  into  K  shows  that  the  (.Cj)  -contribution  to  UN  is  given 
exactly  by  (69)  .  It  should  be  noted  that  this  formula  is 
also  immediately  valid  for  M=1  without  even  introducing  the 
spherical  change  of  coordinates  (73)  ,  since  the  volume  of  the 

unit  interval  is  2  (correctly  given  by  (80)  . 


Although  we  began  the  discussion  by  considering  the  contri¬ 
bution  of  a  single  term  to  the  expectation,  the 


homogeneity  of  the  expression  (74)  shows  that  the  other  terms 
qN  .  2 «  2  *N,  M,2  .  _  ,, 

^1^1  ln  d/l)“entry  will  each  make  an 

identical  contribution.  The  same  argument  applies  to  other 
diagonal  entries  of  and  homogeneity  again  shows  that 

each  term  will  make  the  same  contribution,  given  by  (69)  .  This 
justifies  the  conclusion  that 


-N  mnNi 


(81) 


where  I  is  the  identity  matrix. 

Next,  let  us  transform  back  from  S -coordinates  to  the 
original v -coordinates  as  given  in  (58).  We  have  proven  that 
oux  original  expectations  U^,U2  scalar  multiples  of  the 
covariance  matrix  S: 


uN  =  s"2  (mnNi)  s'2  =  mnN§ 


(82) 


Finally,  substitution  of  (70) -(72)  and  (45)  into 
(69)  followed  by  numerous  cancellations  yields  the  expressions 
(24)-(27)  in  Section  3. 
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APPENDIX  E 


VALIDATION  OF  THE  MODIFIED  RICCATI  EQUATION 


A  set  of  Monte  Carlo  simulations  was  performed  to  validate 
the  result  of  the  modified  (deterministic)  approximation  (33) 
of  the  stochastic  Riccati  equation  (16) . 


The  following  notations  are  used: 


P(klk) 

E|p(k|k)j.- 
(k  Ik)  - 

^(kjk}- 
£d(klk)  - 


The  PDAF-calculated  conditional  error 
covariance  from  the  stochastic  Riccati  equation 
(a  random  matrix,  different  from  run  to  run) 

The  Monte  Carlo  average  of  the  above 

I 

The  true  covariance  of  the  PDAF  estimation 
error  (different  from  run  to  run) 

The  Monte  Carlo  average  of  the  above 

The  result  of  the  modified  (deterministic) 
approximation  (33)  of  the  stochastic  Riccati 
equation 


The  results  are  for  the  numerical  problem  specified  in  [8j. 
where  targets  were  tracked  in  two-dimensional  geographic  space. 
The  sample  means  were  computed  from  ten  independent  runs.  The 
scalar  performance  measures  used  in  comparing  the  matrices  P^ 
e{p}  and  E-fpt}  were 

R.M.S.  error  =  J  tr (£)  (83) 

R.M.S.  position  error  =J  (84) 


Figures  11  -  18  present  the  error  comparison  using  the 


above  two  scalar  measures  for  up  to  200  time  steps  and  a  wide 
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range  of  PD,  Pp.  The  actual  position  error  was  closely 
approximated  by  the  modified  Riccati  equation  even  for  the 
very  low  PD  =  .3.  The  overall  error  was  still  reasonably  well 
approximated  down  to  PQ  -  .4. 

It  appears,  therefore,  that  the  modified  Riccati  equation 
can  be  considered  as  a  reliable  design  tool  except  for  the 
situation  of  very  low  Pp. 
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M=FMS  position  error 


ElPti _ 


Figure  11.  M*Pk|k,  and  M*E{^|k}  for  (PD,Pp)  =  (1,0) 
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: 

: 

* 

3 


U 

o 


Pd  _ 

ElPtJ _ _ 

€CP) _ 


Figure  12-'  M*E{p£|k}f  and  M^EtE^}  for  (PD,Pp)  =  (.9,. 01) 
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M=RMS  error 


Figure  13.  M*fg|k,  M*E{Ejr(k},  and  M*E{£k)k}  for  (PD,Pp)  =  (.8, .02) 
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3 

: 


Figure  14. 


M‘£k|k'  M‘EtEk|k>'  andM-EtP^}  for  (PD,PF)  =  (*7,.03) 
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M-RMS  position  error 


E(PtJ. 


Figure  15.  M*E{p£)k},  and  M’ECE^}  for  (PD,Pp)  =  (.6,. 04) 
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Figure  16.  M-£^|k,  M'E{l£||t},  ^  M‘E{Ek|k}  for  (PD,Pp)  =  (.5,. 05) 
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M-RMS  position  error 


Vi 


p«  _ 


ECPtl.. 


ElPl _ 


Figure  17. 


M‘E{^|k}f  andM-Et^^}  for 


(PD,PF)  =  (.4 ,.06) 
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M=RMS  position  error 


* 

* 


Pd  _ 

ElPtJ _ . 

EIP) _ 


Figure  18.  and  M^EfE^}  for  (PD/Pp)  =  (.3f.07) 

(note  that  scale  of  RMS  error  is  different  from  other  figures) 
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FIGURE  3.  SKETCH  OF  COLSTRIP  UNIT  2  ID  FANS,  DUCTS,  AND  STACK  SHOWING  MEASUREMENT  LOCATIONS. 
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UNIT  2  ARRANGEMENT  IS  MIRROR  IMAGE. 


FIGURE  1.  FAN  MANUFACTURERS  PERFORMANCE  CURVES  FOR  COLSTRIP  UNIT  1  AND  2  ID  FANS  OPERATING  AT  SPECIFIED 
GAS  CONDITIONS. 
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